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1. Summary 


„ . c^ar Magnetic Elements” was a re- 

•h* 

— : - - — th ; ^ 

The primary goa i . n to t he associated magnetic elemen s. magnetic flux elements in 

bright points” in re a 10 ^ bright signatures associate wi * radical in 

the order of 100 km or less in Aw**) ^ A mo lecular bandhead of to ^ ^ 

the photosphere. They «««» resolution and highest con r. 

the solar spectrum and offer the M 0 

magnetic structure on the Sun. , ^ ground -based images from both the 

Although G-band bright points had been P revl °^ ^ ^3 observa tions of Prof. Goran c larmer 
Pic du Midi and Big Bear solar observatone , the attention of many rose 

I the Swedish Vacuum Solar Telescope ( SVST \ ^ br * d brigUt points presented a 

Studies based on G-ba„d brtgh. points ficld . 

to the formation and evolution o , hcsis wor lt of Dr. Berger which estab- 

This investigation was proposed as a follow-on to \ G . band br ight points using the SVS » 

lished most of the ^^^Q^j^i^^deiHnltlarch of 1093 , the three year investigation comprised 
the primary instrument. ungm<u . 
two major programs: 

. ..baerv^^ 

spectra of G-band bright points along ot both magnetic elements (idcnti te 

movies in order to " " ££ s'een G-band. In addition, we study the 
in the magnetogtams) ° c n K-line imaging. 

chromospheric atmospheric .esponse using ^ ^ ^ c)mlctcr , sties 

. A theory /modeling program whose goal is “ sl ™ ^ spoctrum . A full understanding of 
of magnetic elements in the G-band teg, accurate simulation using existing 

Drs. Thomas Berger and Mats Loidahl P « ^ 

currently at the Swedish Solar Ofcsenatore. S‘«^oM Dr L6f(la hl was m charge 

■ phMe 
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. tho SVST Dr. Don Kiselraan of 

arv to achieve diffraction-limited paging at t ical and modeling efforts. 

and observational programs, respecti y observational and theoretical pro- 

‘ Reports " d w,th this " 
o l 

oKc^rvational Milestones: 


• Observing Campaigns at the S 

1. 12 May - 01 June 199S 

2. 24 Aug - 30 Aug 1998 

3. 15 May - 06 Jun 1999 

4. 09 Sep - 16 Sep 1999 

5. 16 May - 25 May 2000 




PPDS Datasets: 

1 ll-Jun-1997: 4 hours o( G-band and 4386 
2 . 12 -Jun- 1997 : 6 hours of G-band and 4386 


A continuum 
A continuum 


‘ G-band,CanK-U„e, and SOUP fiUcrgrams in coord, nation with TRACE, 

and MDI -TRACE MDI coordinated campaign. 

3. 20-May-2000: G-band. LPSP, TRACb, mu 

In terms of scientific output, the W 

the dataset that first elucidated the small seal t0Rr ams in a detailed comparison 

Phenomenon. Using the SOUP filter images, we were able to identify for the first 

with TRACE Fe IX/X 171 k moves and yohkoh/ ^ ^ t0 be the conducttvely 

time the highly dynamic interact, on o bnght mos ^ ^ givt ,„ oss its ret, elated 

TVT A C\X7_Qftnn8 
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dataset was used to produce . r det^d in “ncer. with the 
and the underlying magnetic held using g 

G-band filtergram series. d i 2 June 1997 PPDS datasets represent 

In terms of high resolution photospheric imaging • ^ solar photosp here. Both movies 
the longest and highest resolution movies P I0xim ately 2 arcminute field-of-v,ew 

— - 

TVionretical/Modeling Milestones: 


. Synthesized solar G-band spectra using the SCAN-CH line list 

.1 0. w w» n~ 

,, • the G ban d calculated for 200 km flux tube models. 

• LTE radiation profiles in the G 

. Stokes vector modeling of the SOUP magnetogram response. 

In summary, the observational and ^ magnetized solar 

G-band bright points are a natural y p of this contract lias conclusively dcmon- 

photospheric atmosphere. The obser ™ distin g uishc d from granular edge brightcnings. occur 

strated that G-band bright points, pr P - co-morphous with the magnetized 

only in magnetised regions of the sola, of J 0 , Hydrostatic LTE dux 

structures on Off 5 spatial scales and temporal sc ^ rcptoducc mosl 0 f the relevant contrast 
tube models based on standard atmosp e ^ on6oi „ s in tl ,e area, of non-LTE 

ranges of G- band bright points. However, work b> to determ j ne if there are demonstrable 

radiative transfer and non-standard atmosp answer at the beginning of the 

— While we pointed on, areas that retire further analysis 

r^^derr/the formation and dynam.es o, small-scale magnet, c element, ^ ^ 

A major unexpected contribution of the contract war . the for the study 
understanding of the TRACE mess phenomenon, opening P 
„r dm colnr transition region. 


2. Publications 

w, list here the major publications generated under funding by the contract. The list is sorted m 
chronological order. . 
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, “Small-Scale Magnetic Flux la .he Photosphere: G-band. Marram, and Ca II K-l.ne 

Comparisons” 

M. G. Lofdahl, T. E. Berger, A. Title 

Oral Presentation, AGU Spring Meeting, Boston, 1998. 

2. “Diffusion Measurement Processes 

T. E. Berger, A. M. Title 

Poster Presentation, AGU Spring Meeting, Boston, 1998. 

3, “Joint TRACE and La Palma Observations” 

T. E. Berger, B. de Pontieu, G. B. Scharmer 
Invited talk, AGU Fall Meeting, San Francisco, 1998. 

4 “High Resolution Imaging of the Solar Chromosphere/Corona Transition Region" 

T E Berger, B. De Pontieu, C. J. Schnjver, and A. M. Ti 
Astrophysical Journal Letters, volume 519, pp. L9MOO, 1999. 

5. “Dynamics of Transition Region Moss” 

T E Berger, B. De Pontieu, C. J. Schrijver, and A. M. Title 

Oral Presentation 79.01, AAS Centennial Meeting, Chicago, 2-June-19 • 

6. “Structure and Dynamics of Transition Region Moss" 

PoslrT^ DiaRn ° StiCS m thC SOUr 

Transition Region and Corona 
22-25 June 1999, Paris, France. 

7 ' ;Tb“ De Pontieu. L. Fletcher, T. D. Tarbell, C. J. Schrijvcr. A. NL Title 
Post Plantation, 2nd TRACE Workshop, 24-26 August 1999. Monterey, California. 

8 “What is Nloss? \ \r 'Tif In 

' T. E. Berger, B. De Pontieu.L. Fletcher, T. D. Tarbell. C. J. Schnjver, A. M. Title 
Solar Physics, 190, 409, 1999. 

9. “On the Nature of “Moss” Observed by TRACE” 

P. C. Martens, C. C. Kankelborg, and T. E. Berger 
Astrophysical Journal, 537. 471, 2000. 

10. “Small-scale Magnetic Elements” 

T Berger 

Talk given at the Second Solar-B Science Meeting, ISAS, Japan, Oct. _000. 

■ 11 “On the Relation of G-band Bright Points to the Magnetic Field” . 

T E . Berger and A. M. Title, Astrophysical Journal, accepted for May 10, 2001 issue. 


NASW-98008 



NASW-9800S 


Final Report: Spring 2001 


,, “A Five-hour Sequence of High-Resolution Solar Photospheric Images in Two Wavelengths 

Restored by the use of Phase Diversity 

M.G. Lofdahl, T. E. Berger, and J. H. Seldin, in prepar 

13 . -The formation of G-band bright points I. Standard LTE analysis' D. Kiselman, R. Rntten, 
& Plez, in preparation. 

14. “Spatio-Temporal Multiscaling and Turbulent Solar Flows" 

J. K . Lawrence, A. C. Cadavid, A. A. Rurma.km, and T. E. Berg 
Physical Review Letters, in revision. 

,U addition, the creation of a website detailing the UchheedMarun ground-based so.ar physics 
effort was partially funded by the contract: www.lmsal.com/LaPalma. 


3. Analysis Tools 

Be^r^ 

These codes are written in the ANA language. 




2D Discrete Wavelet Transform (DWT) package. 

Based on the Numerical Recipes in C algorithm, an ANSI-C code was 

using the IDL DLM method. 


developed 


and linked 




Multi-scale image alignment tools. 

IDL code for performing a general affine transformation on images m on «.r to 
alignment and scaling was developed. 


facilitate rigid 




Image Destretching tools. 

Based on the ANSI-C codes of Louis 
code for performing image-to-image 
was developed. 


Strous which are linked to IDL via custom DLMs, IDL 
uniform grid warping as well as time series destretclnng 


Fourier domain filtering tools. , . . svqt 

The contract partially funded the ongoing effort to develop efficient FFT algorithms in ANS - 
C for linking into IDL via the DLM method. In addition. IDL code for performing p-mo e 
supression in the Fourier domain was developed. 


We include the relevant codes as appendices to this report. 
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4. Appendix: Analysis Tool Codes 
4.1. IDL Discrete Wavelet Transform System Routine 

•include <stdio.h> 

•include "export. h" 

•include "nr.h" 

•include "nrutil.h" 

•define NRANSI 

•define ARRLEN(arr) (sizeof (arr)/sizeof (arr [0] ) ) 

•if DLM 

IDL.VPTR dvt(int , IDL.VPTR *); 
int IDL.Load(void) 

static IDL.SYSFUN.DEF function.addr [] = { 

{ dut, "DVT", 3, 3. 0} 

>: 

return IDL.AddSystemRoutineCfunction.addr, IDL.TRUE, 1); 

> 

flendif 

void pwtC float • , unsigned long, int ); 
void print. vptr( IDL.VPTR, char • )l 
typedef struct { 

int ncof , ioff , jof f ; 
float *cc , *cr ; 

} wavefilt; 

uavef ilt wf ilt ; 


IDL.VPTR dwt ( int argc, IDL.VPTR *argv) 

/* Discrete Wavelet Transform. Required inputs: 
a[] : array to be transformed 

filter[]: array of filter coefficients (scaling function), 
isign: isign > 0 does forvard transform, <0 does inverse. 

Returns ndim-dimensional discrete wavelet transform of array a. 
If isign = 1, performs forward transform; if isign = -1 performs 
inverse . 


From Numerical Recipes in C, Press et al . , 2nd edition. With 
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modifications to adapt to XDL: p-t is the onlycallpossible 
(no separate coding of daub4.c or similar functions). IDL ^ 
supplies filter coefficients directly via filtert] - there » 
need for pvtset. NOTE: filter should be unity offset <ie first 
element should be 0) due to use of unity offset Numerical Recipes 

algorithm. 

Wranner routine WAVELET. PRO checks inputs and calls this function. 


********** 




I************************** 


IDL.LONG il, i2, i3, k, n, ncoef , *nn, nnew, ntot=l 
IDL.VPTR filter, rfilter, result; 

int i, scale.iteration, x.dia = 1, y.dim = 2, idim; 
float *a, »b , sig=-1.0, »»ksp; 
short isign; 

UCHAR atype, ndim, rtype; 
extern wavefilt wfilt; 


static IDL EZ ARC arg.structD = { /‘parameter definitions •/ 

{IDL EZ DIM ARRAY, I DL.TYP.B .SIMPLE, IDL.EZ.ACCESS.R, IDL.TYP.UNDEF , 0, 
{IDL~Ez”dIM~ ARRAY , IDL.TYP.B.SIMPLE, IDL.EZ.ACCESS.R, IDL.TYP.UNDEF. 0. 
{IDL.EZ.DIH.HASK(O) . IDL.TYP.B.SIMPLE, IDL.EZ.ACCESS.R, IDL.TYP.UNDEF. 


0) ( /*input array •> 
0), /*filt«r cooffa 
0 # 0), / • i a i gn flag • / 


}; 

IDL_EzCallCleanup(argc, argv, arg. struct); 
IDL.EzCallCargc, argv, arg.struct) ; 


/*get data array •/ 

/•print_vptr(argv[0] , M argv[0] N ) ; •/ 
a = (float *) argv [0] ->value . arr->data ; 
nn = argv [0] ->value . arr->dia , 
ndim = argv [0] ->value . arr->n_ dim ; 


/•get temporary array for result •/ 

b = (float *) IDL_MakeTempArray(IDL_TYP_FLO AT, ndim, nn.IDL_BARR.INI. NOP, Aresult) 


a 

b 

nn 


1 

1 

1 


/•unity offsets from here on out! */ 


/•get filter and load coefficients into external structure (replaces pvtset)*/ 
/•NOTE: filter should have first element - 0 to handle unity offset. Thus 
a Daubechies 4-filter has 5 coefficients: 0,0.482... •/ 

/*print_vptr(argv[l] , M argv [l] ") ; */ 
vfilt.ee = (float *) argv [l] -> value . arr->data ; 

ncoef = (argv [1] ->value . arr->n_elts) -1 ; /*handle unity offset with the -!•/ 
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wfilt.ncof - ncoef ; 

wfilt.cr = (float *) . 

IDL_MakeTempArray ( IDL.TYP .FLOAT , 1 , argv [1] ->valuo . arr->dim, IDL_BARR.INI.ZERO , trf liter ) ; 

for (k=l ;k<=ncoef ;k++) { 

wf ilt . cr [ncoef +l-k] = sig*wf ilt .cc[k] ; 
sig = -sig; 

> 

wfilt.ioff = wfilt.joff = -(ncoef >> 1); 

/♦for (i=0; i<= r ncoef ; i++) printf ('7,f ",wf ilt .cc[i] ) ; 
printf ("\n") ; 

for (i=0;i<=ncoef ;i++) printf ( M 7.f " ,wf ilt .cr [i] ) ; */ 

/* get isign */ 

isign = argv [2] ->value . i ; 

/•printf ( M i3ign=y.d\n'\ isign) ; •/ 

/•Here starts the NR code, (almost) unmodified from original •/ 

/• CHANGES: significant changes have been made so dwt would do wavelet transforms 

• on 2D images correctly. The following code is commented. 

•/ 

for (idim=l ; idim<=ndim; idim++) ntot *= nntidio] ; 
wksp = vector(l.ntot) ; 

for (i=l;i<=ntot;i++) b[i]=a[i]; /*copy input array to avoid overwriting •/ 

/• The purpose of scale.iteration is that at each scale, the effective size of the 

• image is decreased by half, basically quartering each time. In the case of 

• the forward transform, scale iteration is doubled for each iteration, while 

• in the reverse transform, scale iteration starts at a maximum and halves at 

• each iteration. 

*/ 

/•printf ("This is the beginning . \n" ); •/ 

if (isign >= C) { /•forward wavelet transform condition •/ 
scale.iteration = 1; 

} else { /*reverse wavelet transform condition*/ 
scale.iteration = nn[x_dim]/8; 

} 

/•printf ("hello\n initial scale: */,i\n M , scale.iteration) ; */ 

while (scale.iteration > 0 fet nn [x_dim] /scale.iteration >= 8 kk nn [y.dim] /scale. iteration >= 8 
/•printf ("scale.iteration: 7,i\n" , scale. iterat ion) ; •/ 

n = nn [x_dim] /scale.iteration ; 
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for (i2 = 0; i2 < ntot/scale.iteration; i2+=nn[y_dim] ) { 
for (i3 = i2+l, k*l; k <= n; k++, i3 += 1) { 
wksp[k] = b[i3] ; 

> 

pwt(wksp,n,isign); 

for (i3 = 12+1, k=l; k <= n; k++, i3 += 1) { 
b[i3]=wksp[k] ; /‘copy back from workspace •/ 

> 

> 

/•printf ("part l\n");*/ 
n = nn[y .dim] /scale.iteration; 

for (il = 1; il <= nn[y_dim] /scale.iteration; il += 1) { 
for (i3 = il, k=l; k <= n; k++, i3 += nn[x_dim]) { 
wk3p[k] = b[i3] ; 

> 

pwt(wksp,n,isign) ; 

for (i3=il, k=l; k <= n; k++, i3 += nn[x_dim]) { 
b[i3]=wksp[k] ; /»copy back from workspace • / 

> 

> 

/•printf ("part 2\n");*/ 

if (isign >= 0) { /‘forward wavelet transform condition • / 
scale.iteration •= 2; 

> else { /‘reverse wavelet transform condition*/ 
scale.iteration /= 2; 

> 

/•printf ("scale.iteration at the end is : 7.i\n", scale iteration) • •/ 

> 

/•printf ("good\n") ;•/ 
free_vector(wksp, 1 ,ntot) ; 

IDL_EzCallCleanup(axgc , argv, arg.struct) ; 

/•printf ("even better\n") ; •/ 
return result; 


void pwt(float a[] , unsigned long n, int i 3 ign) 

/• Partial wavelet transform. Applies an arbitrary wavelet 
filter to data vector a[l...n] passed in a[] for isign=l; 
for isign=-l, applies inverse transform. 

From Numerical Recipes in C, 2nd edition. 


{ 


/ 
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float ai ,ail ,*wksp; 

unsigned long i,ii, j , jf i jr t k,nl,ni p nj ,nh,nmod; 
extern wavefilt wfilt; 

if (n < 4) return; 
wksp=vector (l,n) ; 

nmod=wf ilt .ncof *n; /* A positive constant equal to zero mod n */ 
nl^n- 1 ; /* Mask of all bits, since n is a power of 2 •/ 

nh=n >> 1; 

for (j=l; j<=n; j++) wksp[j]=0.0; 
if (isign >= 0) { /* Apply filter */ 

for (ii=l,i=l;i<=n;i+=2,ii++) { 

ni=i+nmod+wf ilt . iof f ; /•Pointer to be incremented and wrapped around •/ 

nj=i+nmod+wf ilt * joff ; 

for (k=l ;k<=wf ilt .ncof ;k++) { 

jf=nl It (ni+k); /*we use bit-wise AND to wrap around the pointers •/ 

jr=nl It (nj+k) ; 

wksp[ii] += wf ilt .cc[k] *a[jf +1] ; 
wksp[ii+nh] += wf ilt .cr [k] •atjr+l] ; 

> 

} 

} else { /* Apply transpose filter •/ 

for (ii=l,i=l;i<=n;i+=2,ii++) { 
ai=a[ii] ; 
ail=a[ii+nh] ; 
ni=i+nmod+vf ilt . iof f ; 
nj=i+mnod+wf ilt . jof f ; 
for (k=l ;k<=vf ilt .ncof ;k++) { 
jf=(nl k (ni+k))+l; 
jr=(nl k (nj+k))+l; 
wksptjf] += vf ilt . cc [k] *ai ; 
wksp[jr] += wf ilt .cr [k] *ail ; 

> 

} 

} 

for ( j=l ; j<=n; j++) a[j]=wksp[j] ; /*copy the results back from workspace •/ 
f ree_vector(wksp, 1 ,n) ; 

} 

#undef NRANSI 


FUNCTION wavelet, array, filter, isign, OVERWRITE=overwr ite 
; Wrapper function for wavelet transform implemented in C. 
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;Applies the wavelet filter FILTER, specified by a string 
j corresponding to values in WT_FILTER.PRO, to the N~dimensional 
; array passed in ARRAY. 

; DVT is composed of c-code routines vtn and pvt adapted, 

; and modified heavily, from Numerical Recipes 2nd edition. 


0N_ERR0R,2 

if N _P ARAMS () It 3 then begin 

MESSAGE, 'Wavelet transform requires three arguments: array, filter string, and isign* 
RETURN,-! 
end 

; Check ARRAY: 
sz = SIZE(axray) 

badnevs = 'All dimensions of ARRAY must be povers of 2* 
for i=l,sz(0) do $ 
if not ispov2(sz(i)) then begin 
MESSAGE.badnevs 
RETURN, -1 

end 

if sz(3) ne 4 then array = FLOAT(array) ;This step is crucial since the C-code uses 
; sizeof (float) offsets to do the transform steps. 

;Check filter: Numerical Recipes algorithm needs unity offset vectors: 
filter = STRUPCASE (filter) 
coef f =vt .filter (filter) 

if coef f (0) ne -1 then filt = [O.,coeff] else RETURN , coef f 
nn = [sz ( 1 ) ] 

for i=2,sz(0) do nn = [nn,sz(i)] 
nmin = MIN (nn) 

if N_ELEMENTS(f ilt)-l gt nmin then begin 

MESSAGE , * Filter too long for data array in at least one dimension’ 

RETURN , - 1 
end 

;help,f ilt 

;for i=0,N_ELEMENTS(f ilt)-l do print , f ilt (i ) 

;Check isign 
case 1 of 

isign gt 0: isign = 1 
isign It 0: isign = -1 
else: begin 
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MESSAGE, 1 specify ISIGN > 0 for forward transform, ISIGN < 0 for inverse 1 
RETURN, -1 
end 
end 

;A11 okay - do transform in C: 

result = dwt (array, FLOAT(filt) ,isign) ;THE FLOAT(filt) IS CRUCIAL!! 

if KEYWORD. SET(overwrite) then begin 
array = result 
RETURN, 1 

end else RETURN, result 
END 


FUNCTION wt.filter, filter 


; Look-up table cased on string FILTER. Returns the 
; wavelet filter coefficients in a double array COEFFS. 

;Coeff icients gotten from Mathematica Wavelet Explorer (varying precision ones) 
; and from Numerical Recipes (const, precision ones). 


filter = STRUPCASE (filter) 

CASE filter of 

’HAAR’: coeff = R£PLICATE( 1 . /SQRT(2 .0D) ,2) 

'DAUB4 * : coeff « FL0AT( [0.4829629131445341 ,$ 

0.8365163037378079,$ 

0.2241438680420134,$ 

-0.1294095225512604]) 

’DAUB6': coeff = FL0AT( [0.3326705529500828 ,$ 

0.8068915093110929,$ 

0.4598775021184916,$ 

-0.1350110200102547,$ 

-0.08544127388202671,$ 

0.03522629188570956]) 

’DAUBS’: coeff = FL0AT( [0.2303778133088966 , $ 

0.7148465705529158,$ 
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0.6308807679298591,$ 
-0.02798376941686012,$ 
-0.1870348117190932,$ 
0.03084138183556081,$ 
0.03288301166688522,$ 
-0.01059740178506904] ) 

’DAUB10’ : coeff = FL0AT([0. 1601023979741925,$ 

0.6038292697971881,$ 
0.7243085284377716,$ 
0.1384281459013218,$ 
-0.2422948870663802,$ 
-0.03224486958463779,$ 
0.07757149384004565,$ 
-0.006241490212798174,$ 
-0.01258075199908195,$ 
0.003335725285473758] ) 

’DAUB12’ : coeff = FLOAT([0. 111540743350,$ 

0.494623890398,$ 

0.751133908021,$ 

0.315250351709,$ 

-0.226264693965,$ 

-0.129766867567,$ 

0.097501605587,$ 

0.027522865530,$ 

-0.031582039318,$ 

0.000553842201,$ 

0.004777257511,$ 

-0.001077301085]) 

’DAUB20’ : coeff = FL0AT( [0.026670057901 ,$ 

0.188176800078,$ 

0.527201188932,$ 

0.688459039454.$ 

0.281172343661,$ 

-0.249846424327,$ 

-0.195946274377,$ 

0.127369340336,$ 

0.093057364604.$ 

-0.071394147166,$ 

-0.029457536822,$ 

0.033212674059,$ 

0.003606553567,$ 

-0.010733175483,$ 

0.001395351747,$ 
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0.001992405295,$ 

-0.000685856695,$ 

-0.000116466855,$ 

0.000093588670,$ 

-0.000013264203]) 

’MEYERO' : coeff = FL0AT( [0.003691882686403731 ,$ 

0.004541394739570804,$ 

-0.02000702924793551,$ 

0.0225516658829227,$ 

0.04050242501133709,$ 

-0.1000351462396785,$ 

-0.057691904077159.$ 

0.4316322262631767,$ 

0.771509959510067,$ 

0.4316322262631767,$ 

-0.057691904077159,$ 

-0.1000351462396785,$ 

0.04050242501133709,$ 

0.0225516658829227,$ 

-0 . 0200070292479355 1 , $ 
0.004541394739570804,$ 
0.003691882686403731]) 

’ MEYER 1 ’ : coeff = FL0AT([O. 01260563175115622,$ 

-0.01472281058330518,$ 

-0.02431045382096553,$ 

0.04603894635662135,$ 

0.03661850739611902,$ 

-0.1194679203378185,$ 

-0.04603454889385984,$ 

0.4391769769258334.$ 

0.7566719914327569,$ 

0.4391769769258334.$ 

-0.04603454889385984,$ 

- 0 . 1194679203378185 ,$ 

0.03661850739611902,$ 

0 . 04603894635662135 ,$ 

- 0 . 02431045382096553 ,$ 

- 0 . 01472281058330518 ,$ 

0 . 01260563175115622 ]) 


else: begin 

PRINT, ’Unsupported filter type in WT.FILTER.PRO * 
PRINT, ’Supported filter types:’ 

PRINT,’ Haar’ 
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PRINT,' Daub4 , Daub6, Daub8, DaublO, Daubl2, Daub20’ 
PRINT,* MeyerO, Meyer 1' 
coeff = -1 
end 


END 

RETURN, coeff 
END 


4.2. IDL Image Alignment Tools 


FUNCTION AFFINE, image, ax, my, ax, theta, xc, yc,$ 

INTERP=interp, CUBIC=cubic, MISSINC=mi3sing 


NAME: 

AFFINE 

PURPOSE: 

Apply the affine transf ormation given by the input parameters 
to IMAGE. 


CATEGORY: 

Z3 - Image processing, geometric transforms, image registration. 

CALLING SEQUENCE: 

transf ormed_image = AFFINE (image ,mx , my , sx , theta , xc , yc) 

INPUTS: 

IMAGE: The image to be transformed. Must be 2-D. 

MX, MY: Magnification factors in x and y directions. 

SX: Horizontal shear term. 

THETA: Rotation angle in DEGREES. THETA > 0 => counterclockwise rotation. 

XC, YC: Center of rotation if THETA ne 0 else translation. 

KEYWORDS : 

INTERP: Set this keyword for bilinear interpolation. If this keyword 
is set to 0 or omitted, nearest neighbor sampling is used. 
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Note that setting this keyword is the same as using the 
ROT.INT User Library function. This change (and others) 
essentially makes ROT.INT obsolete. 

CUBIC: If specified and non-zero, "Cubic convolution" 

interpolation is used. This is a more 

accurate, but more time-consuming, form of interpolation. 

CUBIC has no effect when used with 3 dimensional arrays. 

If this parameter is negative and non-zero, it specifies the 
value of the cubic interpolation parameter as described 
in the INTERPOLATE function. Valid ranges are -1 <= Cubic < 0. 
Positive non-zero values of CUBIC (e.g. specifying /CUBIC) 
produce the default value of the interpolation parameter 
which is -1.0. 

MISSING: The data value to substitute for pixels in the output image 
that map outside the input image. 

OUTPUTS: 

NONE 

RETURNS: 

TIMAGE: the affine transformation of input image IMAGE. 

COMMON BLOCKS: 

None . 

SIDE EFFECTS: 

None . 

RESTRICTIONS: 

None . 

PROCEDURE: 

Uses POLY. 2D to warp the input image according to the 
given parameters. 

See: Image Processing for Scientific Applications 

Bernd JV'ahne 

CRC Press, 1997, Chapter 8. 

Same as ROT. PRO but includes shear term and /PIVOT is assumed. 

MODIFICATION HISTORY: 

T. Berger, LMATC, 26-February-1998 . 

T. Berger, LMATC, 25-June-1998: made THETAR internal variable. 
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0N.ERR0R.2 

sz = SIZE(image) 
if sz[0] ne 2 then begin 

MESSAGE, * Input image must be 2-D* 

RETURN , - 1 

end 

if N.PARAMSO eq 1 then mx = 1.D0 else mx = DOUBLE(mx) 

if NLPARAMSO le 2 then my = 1.D0 else my = DOUBLE(my) 

if N.PARAMSO le 3 then sx = 0.D0 else sx = DOUBLE(sx) 

trans=0 

if N_PARAMS() le 4 then begin 
theta = 0.D0 
trans = 1 

end else thetar = DOUBLE(theta* ! DTOR) 

if thetar eq 0.0 then trans=l 
if N.PARAMSO le 5 then begin 
if thetar eq 0 then begin 
xc = 0 
yc = 0 

end else begin 

xc = D0UBLE(sz[l]/2.) 
yc = DOUBLE (sz [2] /2.) 

end 

end else xc = DOUBLE(xc) 
if N.PARAMSO le 6 and thetar eq 0 then yc 
if N_PARAMS() le 6 and thetar ne 0 then yc 


= 0 

= D0UBLE(sz[l]/2. ) else yc = DOUBLE(yc) 


axy = mx*my 


if trains then POO = xc else $ 

POO = ( (-my*xc + mx.sx*yc).COS(thetar) + mxyxc + mx.yc.SIN(thetar) )/=uy 
P10 = -sx*C0S (thetar) /my - SIN(thetar) /my 
P01 = COS(thetar)/mx 
Pll = 0. 

P = [P00,P10,P01,P11] 


if trans then Q00 = yc else $ 

Q 00 = ( mxy.yc - mx.ycCOS(thetar) + (-my.xc*mx«sx*yc) .SIM(th«tar ) )/mxy 
Q10 - COS (thetar) /my - sx*SIN(thetar) /my 
COl = SIN (thetar) /mx 
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Qll = 0. 

q = [qoo.Qio.Qoi.qii] 

i=0 

if KEYWORD_SET(interp) then i=l ;bilinear 
if N_ELEMENTS (cubic) eq 0 then cubic = 0 

if N.ELEMENTS (missing) eq 0 then return, P0LY_2D( image, P, Q, i ,CUBIC=cubic) else $ 
return , P0LY.2D ( image , P , Q , i ,CUBIC=cubic , MISSING=mi ssing) 

END 


PRO AFFINE_SOLVE, xin, xpin, $ 

mx, my, sx, theta, xc, yc,$ 
VERBOSE=blab 


NAME: 

AFFINE.SOLVE 

PURPOSE: 

Calculate the parameters of a general affine image 
transformation given a set of points from tuo images: 
one of the images is assumed to be the reference image, 
the other is assumed to be an image translated, rotated, 
scaled, and possibly sheared relative to the reference image. 

The form of the general transformation is affine: 

X = tranformed coordinates = [T+ HSR T-] X' 
where, in homogeneous coordinates, 

* = TRANSPOSE [x, y, 1] ; test image vector 

T+ = [[1,0, xO] ,(0,1, yO] , [0,0,1]] : translation of (0,0) back to (iO yO) 
M = [[mx ,0,0] , [0,my ,0] , [0,0, 1]) : scale 
S - [[1 , sx ,0) , [0 , 1 ,0] , [0,0 , 1] ] : horizontal shear 

R = [[cos(t) ,-sin(t) ,0] , [sin(t) ,cos(t) ,0] , [0,0, 1]] : 

rotate clockwise by angle t about origin. 

T - [ [1 ,0 , -xO] , [0, 1 ,-y0] , [0,0, 1]] : center of rotation to (0,0) 

X' = TRANSPOSE [x ' , y 1 , 1] : reference image vector 

CATEGORY: 

Imag. processing, geometric transforms, image registration. 
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CALLING SEQUENCE: 

AFFINE_SOLVE , xin , xref in , sx , sy , s , theta , xO , yO 
INPUTS: 

XIN: 2xN dimensional array of points taken from imagel 

which correspond to the same points in the reference image. 
Xi = XIN(0,*) 

Yi = XIN(1,*) 

N is the number of points. 

XPIN: 2xN dimensional array of points from the "reference image" 

which correspond to points in the image. 

KEYWORDS: 

VERBOSE: If set, print the transformation elements to the screen. 
OUTPUTS: 

MX, MY: Magnification factors in x and y axes, respectively. 

SX: Horizontal shearing factor. 

THETA: Rotation angle in degrees. 

XC,YC: Center of rotation vector elements OR translation 

vector elements. 


COMMON BLOCKS: 

None . 

SIDE EFFECTS: 

None. 

RESTRICTIONS: 

N, the number of matched points in the transformed and reference 
images should be large (greater than 20), should be taken from 
widely spaced locations in the image f ield-of -view , and should 
be measured to within 1-pixel for greatest accuracy. 

Off-center rotation and translation require a two-stage approach 
for image registration, i.e. in the first stage, apply the parameters 
given by this routine to the test image. A second set of points 
is then selected from the image and the reference image, and 
a second run of this program should output a final translation 
to be applied to the test image to bring it in registration with 
the reference image. This is tested for and the user is alerted. 
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; PROCEDURE: 

; Using least squares estimation, determine the elements 

; of the general affine transformation (rotation and/or scaling 

; and/or translation and/or shearing) of an image onto a reference 

; image . 

; See: Image Processing for Scientific Applications 

I Bernd JV'ahne 

’> CRC Press, 1997, Chapter 8. 

; Use AFFINE. PRO (or ROT. PRO if no shear is found) to apply the 

; transformation to the test image after computing them with this routine 

; MODIFICATION HISTORY: 

; Written: T. Berger, LMATC, 24-Feb-1998. 

; Added no rotation/translation test. TEB, 2-March-98. 

0N.ERR0R.2 

x = xin 
xp = xpin 


nl = (SIZE(x)) [1] 
n2 = (SIZE(xp)) [1] 
if n2 ne nl then begin 

MESSAGE, 'Number of points must be the same in both input arrays 
RETURN 


end 


y = D0UBLE(x [* , 1] ) 
x = DOUBLE (x[*,0]) 
yp = D0UBLE(xp[* , 1 ] ) 
xp = D0UBLE(xp[*,O] ) 


;Least squares solution for matrix elements: see Notebook 5 p 132 
AT = DBLARR(2*nl ,6) 
zerow = [REPLICATE(0,nl)] 
onerow = [REPLICATE(1 ,nl)] 

AT [INDGEN (nl ) *2 , •) = [[x] , [y] .onerow, zerow.zerow, zerow] 
AT[INDGEN(nl) *2+1 , *] = [zerow, zerow.zerow , [x] , [ y ] , onerow] 

A = TRANSPOSE(AT) 

b = TRANSP0SE(DBLARR(2*nl) ) 

b[I.TOGEN(nl)*2] = xp 

b [INDGEN (nl) *2+1] = yp 

xb = INVERT (AT## A , /DOUBLE) ##ATS#REFORM(b) 
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;Solve for transformation elements: 
theta = -ATAN(xb[3],xb[4]) 

mi = xb [3] * (xb [1] *xb [3] -xb [0] *xb [4] ) /SIN (theta) / (xb [3] *2+xb [4] *2) 
my = -xb [3] /SIN (theta) 

sx = (xb [0] *xb [3] +xb [1] *xb [4] ) /(xb [0] *xb [4] -xb [1] *xb [3] ) 

translation solution: 

det = xb [0] *xb [4] -xb [1] *xb[3] 

denom = (xb[0]+xb[4] ) - det - 1. 

xc = - (xb [2] - xb [2] *xb [4] +xb [1] *xb [5] ) / denom 

yc = -(xb [2] *xb [3]+xb[5] -xb [0] *xb[5] )/denom 

trans=0 

if ABS(theta) It 5e-03 then begin ;no rotation - use simple translation solve: 
trans=l 
xc = xb[2] 
yc = xb[5] 

end 

;Return the transformation FROM x TO xp: 

; ie. rotate image by theta degrees counterclockwise to get to reference image 
theta = theta/ Idtor 

if theta gt 0 then dir=’clockvise’ else dir=’counterclockuise’ 


if KEYWORD.SET (blab) then begin 

stheta = STRCOHPR£SS(ABS (theta) ,/re) 

sxc = STRCOMPRESS (xc , /re) 

syc = STRCOMPRESS (yc, /re) 

smx = STRCOMPRESS (mx, /re) 

smy = STRCOMPRESS (my , /re) 

ssx = STRCOMPRESS (sx, /re) 

PRINT, ’ ’ 

PRINT, ’Relative to the reference image, the test image is:’ 


if not trams then begin 

PRINT,’ Rotated stheta, ' 

PRINT,’ vith the center of 

PRINT,' xc = \sxc 

PRINT,’ yc - * f syc 

end else begin 


degrees ’ ,dir 
rotation at ’ 


PRINT, ’ 
PRINT, ’ 
PRINT, ’ 
PRINT, ’ 

end 


Rotated ’, stheta,' degrees ’ ,dir 
Translated by' 

dx = ' ,sxc 
dy = ’ , syc 
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PRINT,' Scaled by a factor of ’,smx,’ horizontally' 

PRINT,’ Scaled by a factor of ’,smy,’ vertically’ 

PRINT,’ Sheared horizontally by a factor of = ’,ssx 
PRINT, ' ’ 

end 

RETURN 

END 


NAME 

IMAGE-REGISTER 

PURPOSE: 

Widget-based program to determine the geometric registration 
between two images of the same scene. Uses affine transformations 
to calculate relative scaling, rotation, displacement, and shear. 

CALLING SEQUENCE: 

image.register , iml , im2 ,r ,xp 
INPUTS: 

iml - image 1, defined as the "reference image" 

im2 - image 2, image to be aligned/scaled/rotated to image 1. 

OUTPUT: 

x - vector of points [2,N] where [0,N] are the x-values of 
reference tie points selected by the user and [1,N> are 
the y-values. 

xp - vector of identical points located by the user in image2. 
KEYWORDS : 

center - not implemented yet. 

output STRING filename. If set, program writes X and XP to 
the file specified by the string filename. 

USEAGE: 

Upon calling, a widget with four display areas opens. In the 
left half of the widget the reference image (iml) is displayed 
with a 4X zoom in the zoom window. In the right half. im2 is 
similarly displayed. The user selects a point in iml using 
the mouse and the left button. The user then finds this exact 
same point in im2 and selects it using the mouse. For a robust 
solution, usually more than 10 points axe required. 

If an incorrect selection is made, use the EDIT>Undo Selection 
menu item. 
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To see the points selected use the DISPLAY>selections menu 
item. 

To print the current solution in the active IDL window, use 
the TRANSFORM>aff ine menu selection. Repeat this step periodically 
during the selection process to see if the solution is converging. 
In some cases, mis-paired points may cause the solution to 
converge to erroneous values; this requires a restart of the 
procedure . 


NOTES : 

Uses AFFINE. SOLVE. PRO to determine the affine transformation 
REVISION HISTORY: 

Written T. Berger 20-April-1998 


PRO Image. Register. event , event 


The main event handler called by Manager, "event" is a 
structure returned by a the WIDCET. EVENT function within 
XManager . 



COMMON im_register, $ 

iml , iml.xO, iml,yO, $ 

im2, im2_x0, im2_y0, $ 

xyvisl, xyvis2, % 

magnl, magn2, $ 

draw lid , draw2id, $ 

detlid, det2id , $ 

winsz, $ 

outf ile ,unitnum, $ 

mousef lag , xymouse , $ 

xl » yl, xlf lag , $ 

x2, y2 , x2f lag , $ 

xr > xpr, $ 

detsz 


; Image 1 and its dimensions. 

; Image 2 and its dimensions. 

;Coordinates of draw windows (LL corner). 
;Current magnification of images 1 and 2. 
;Vindow IDs for display vidgets. 

; Vindov IDs for detail displays. 

;Size of the draw windovs. 

;Qutput file flag and unit number. 

» Flag and position of mouse. 

,Selected points in image 1 and flag. 
;ditto for image 2. 

; Returned arrays. 

;Size of detail windows 


COMMON widget, ids, $ 

irbase , $ 

drawl, drav2 , $ 

detain, detail2, $ 

text lx , text ly , $ 

text2x, text2y 


;Main base ID. 

; Drawing areas. 

;Detail view areas. 

;Text output areas for image 1. 
;ditto for image 2. 
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VIDCET^CONTROL, event. id, GET_UYALUE=uval 

;print , * Event : event 

; print, ’uval = 1 , uval 

case uval of 

'menu 1 : case event. value of 
’File .Save* : begin 


end 

’File. Save As 1 : begin 


end 

’File. Exit': begin ;close the points pairs output file if still open 

if outfile eq 1 then CLOSE, unitnum 
VIDGET.CONTROL , event. top, /DESTROY 

end 

’Edit .Undo Selection’ : begin 
if xlflag eq 1 then begin 
nel = (SIZE(xl) ) [l] 
badx = xl [nel-1] 
bady = yl [nel-1] 

PRINT, badx, bady , ' : Selection removed. Reselect in Imave 1.' 
xl = xl [0:nel-2] 
yl = yl[0:nel-2] 
x2f lag = 1 
end else begin 

nel = (SIZE(x2) ) [l] 
badx = x2 [nel-1] 
bady = y2[nel-2] 

PRINT, badx , bady , . Selection removed. Reselect in Imave 2 ’ 
x2 = x2[0:nel-2] 
y2 = y2[0:nel-2] 
xlflag = 1 

end 

end 

’Display. Selections. Print’ : begin 
npts = (SIZE(xl) ) [l] 
if npts gt 1 then begin 
xlt = xl [l : *] 

ylt = yl [1 : *] 
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x2t = x2[l:«] 
y2t = y2[l :*] 

PRINT,' ’ 

PRINT, ’Selections so far:’ 

PRINT,” 

PRINT,’ IMAGE 1 IMAGE 2’ 

for i=0,npts-2 do begin 

sxl = STRCOMPRESS(xlt[i] ,/re) 
syl = STRCOMPRESS(ylt [i] ,/re) 
sx2 = STRC0MPRESS(x2t[i] ,/re) 
sy2 = STRC0MPRESS(y2t[i) ,/re) 

PRINT, ’( ’.sxl,’, ’ ,syl, ’ )’,’ ( ’,sx2,’, ’,sy2,’ )> 

end 

PRINT,” 

end else PRINT, 'No points selected yet.' 

end 

'Display. Selections. Circle’ : begin 
nptsl = (SIZE(xl) ) [l] 
npts2 = (SIZE(x2))[l) 
if nptsl gt 1 and npts2 gt 1 then begin 
xlt = xl[l:») 
ylt = yl[l:»] 
x2t = x2[l:*] 
y2t = y2[l :*] 

WSET, drawl id 

for i=0,nptsl-2 do TVCIRCLE, 10, xlt [i] «magnl ,ylt [i] •magnl ,color*200 
WSET.draw2id 

i-0 ,npts2-2 do TVCIRCLE, 10, x2t [i] *nagn2, y2t[i] •nagn2,color m 200 
end else PRINT.’Not enough points selected yet.’ 

end 

’Display. Selections. Highlight’ : begin 
nptsl = (SIZE(xl) ) [1] 
npts2 = (SIZE(x2) ) [l] 

if nptsl gt 1 and npts2 gt 1 then begin 
xlt = xl [1 : *] 
ylt = yl [1 : »] 
x2t = x2 [1 : •] 
y2t = y2 [1 : •] 

WSET .draw lid 

for i=0,nptsl-2 do PLOTS, [xlt [i] .magnl] , [ylt [i] .magnl] , $ 
psym=4 , thi=2 , s 733 = 0 . 2»magnl , color=200 , /dev 
WSET,draw2id 

for i=G r npts2-2 do PLOTS , [x2t [i] *aagn2] , [y2t [i] •magn2] , $ 

psym=4 l thi=2,Sjms=0.2*magn2,color=200 f /dev 
end else PRINT, ; Not enough points selected yet. 1 

end 
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'Display. Color Table. Xloadct ...' : XLOADCT 
'Display. Color Table. XPalette. : XPALLETE 
'Transform. Affine’ : begin 
npts = N.ELEMENTS(xl) 
if npts gt 1 then begin 
xlt = xl [1 : •] 
ylt = yl [1:*] 
x2t = x2 [1 : •] 
y2t = y2[i:*] 
x = [[xlt] , [ylt]] 
xp = [ [x2t] , [y2t]] 
af f ine_solve ,x,xp,/verb 
end else PRINT, 'No points selected yet.' 

end 

end 

'menul': case event. value of 

' Zoom. 1:4': begin 

IR.redraw, drawl, 0.25 

end 

' Zoom. 1:2*: begin 

IR.redraw, drawl, 0.5 

end 

' Zoom. 1:1': begin 

IR.redraw, drawl, 1.0 

end 

' Zoom. 2 : 1 * : begin 

IR.redraw, drawl, 2.0 

end 

'Zoom. 4:1': begin 

IR.redraw, drawl, 4.0 

end 

'Zoom. 8:1': begin 

IR.redraw, drawl, 8,0 

end 

'Move. Image Center ' : begin 

xyvisl [0] = iml_x0/2»aagnl-uinsz/2 
xyvisl [ 1 ] = iml_y0/2«nagnl-vinsz/2 
WIDGET.CONTROL, drawl ,SET_DRAW_VIEW=xyvi 3 l 

end 

'Move. Mouse. . . ' :begin 

PRINT, 'Image centers on next mouse down...' 
mouseflag = l 

end 

end 
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'menu2': case event. value of 

'zoom. 1:4* : begin 

IR_redraw, drav2, 4 

end 

'Zoom. 1 :2' : begin 

IR.redraw, draw2, 0.5 

end 

' Zoom. 1:1*: begin 

IR.redrav, drav2, 1.0 

end 

'Zoom. 2:1': begin 

IR.redraw, draw2, 2.0 

end 

' Zoom. 4 : 1 * : begin 

IR. redraw, draw2, 4.0 

end 

' Zoom. 8 : 1 * : begin 

IR.redraw, draw2 # 8.0 

end 

'Move. Image Center 9 :begin 

xyvis2[0] = im2_xO/2»aagn2-vinsz/2 
xyvis2[l] = im2_y0/2*nagn2-vin3z/2 
WIDGET_C0NTR0L,draw2,SET_DlUW_VIEV=xyvi32 

end 

'Move. Mouse. . . ' :begin 

PRINT t 'Image centers on next mouse down... 1 
mouseflag « 1 

end 

end 

'drawl 1 ; begin 

xm = FLOAT(event .x)/magnl 
ym = FLOAT(event .y)/magnl 

case event. type of 

begin ;only register point if previous pair is done and move. mouse is inactive 
if x2f lag eq 1 and mouseflag eq 0 then begin 
x2f lag = 0 

PRINT, ’Image 1: ’ .rrn.ym 
xl = [xl.xm] 
yl = [yi.ym] 
xlf lag = 1 

end else if mouseflag eq 0 then $ 

PRINT, 'Please select the previous point in Image 2 .’ 
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end 

2: begin 

IR_detail,detaill, xm, ym 
strx = STRCOMPRESS(xm,/re) 
stry = STRCOMPRESS (ym , /re ) 

WIDGET.. CONTROL , text lx , SET_VALUE=st rx 
WIDGET_CONTROL,textly l SET.VALUE=stry 

end 
else : 


end 

end 

, drav2 , : begin 

xm = FLOAT (event .x)/magn2 
ym = FLO AT(e vent .y)/magn2 
case event. type of 

0:begin 

if xlflag eq 1 and aouseflag eq 0 then begin 
xlflag = 0 

PRINT, linage 2: 1 ,xm,ym 

PRINT, * * 

x2 = [x2,xm] 

y2 = [y2,ym] 

x2f lag = 1 

end else if mouseflag eq 0 then $ 

PRINT, 'Please select a nev point in Image 1.' 

end 

2: begin 

IR.detail ,detail2 , rm f ym 
strx = STRCOMPRESS (za, /re) 
stry = STRCOMPRESS(ya,/re) 

VIDGET.CONTROL , text2x , SET.VALUE=strx 
WIDGET.CONTROL , text2y , SET_VALUE=st ry 

end 

else : 
end 

end 


end ;case uval 
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RETURN 

END 

PRO IR.Redraw, vid, mag, XCENTER=xcen, YCENTER=ycen 
; Redispalys image after changes to MAG or ORIGIN 


COMMON im^register, $ 

ial, iml.xO, iml.yO, $ 

im2, im2_x0, im2_y0, $ 

xyvisl, xyvis2, $ 

magnl, magn2, $ 

drawlid,drav2id, $ 

detlid, det2id, $ 

winsz, $ 

outf ile ,unitnum, S 

mousef lag,xymouse, $ 

xl, yl , xlf lag, $ 

x2, y2, x2flag, $ 

xpr, $ 

detsz 


; Image 1 and its dimensions. 

; Image 2 and its dimensions. 

Coordinates of draw windows (LL corner). 
; Current magnification of images 1 and 2. 
;Vindow IDs for display widgets. 

;Vindov IDs for detail displays. 

;Size of the draw windows. 

;0utput file flag and unit number. 

;Flag and position of mouse. 

^Selected points in image 1 and flag. 
;ditto for image 2. 
i Returned arrays. 

;Size of detail windows 


COMMON widget.ids, $ 

irbase, $ 

drawl, drav2, $ 

detaill, detail2, $ 

textlx , textly , $ 

text2x , text2y 


;Main base ID. 

Craving areas. 

;Detail view areas. 

;Text output areas for image 1. 
;ditto for image 2. 


if N_ELEMENTS(mag) eq 0 then mag = 1 
if vid eq draul then begin 
vnum = drawl id 
xyvis = xyvisl 
factor = mag/magnl 
zim = iml 
xs = iml_xO 
ys = iral.yO 
end else begin 

wnum = draw2id 
xyvis = xyvis2 
factor = mag/magn2 
zim = im2 
xs = im2_x0 
ys = im2_y0 
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end 

;find center of window on current image or set request: 

VIDGET ..CONTROL, wid, GET_DRAV..VIEV=xyvis 

if KEYWORD_SET(xcen) then xc = xcen else xc = xyvis [0]+winsz/2 
if KEYWORD_SET(ycen) then yc = ycen else yc = xyvis [1] +winsz/2 
;PRINT,xc,yc 

; display zoomed image 

VIDGET_C0NTR0L , wid, DRAV_XSIZE=xs*mag, DRAV_YSIZE=ys*mag 
VSET, wnum 

VIDGET.CONTROL , /HOURGLASS 
TVSCL, conf ac(zim,mag) 
xc = xc*factor 
yc = yc*factor 

xyvis = [xc-winsz/2,yc-winsz/2] 

VIDCET.CONTROL, wid, SET_DRAV_VIEV=xyvis 

jupdate common block: 
if wid eq drawl then begin 
magnl = mag 
xyvisl * xyvis 
end else begin 
magn2 = mag 
xyvis2 = xyvis 

end 

RETURN 

END 


PRO IR_Detail, wid, xa, ym 

Dispalys magnified image in detail window. 
Constant update from mouse notion in main window. 


COMMON im_register, $ 

iml , iml_xO, iml_yO, $ 

im2, im2_x0, im2_y0, $ 

xyvisl, xyvis2 , $ 

magnl, magn2, $ 

drawlid,draw2id , $ 

detlid, det2id , $ 

winsz, $ 

outf ile ,unitnum, $ 

mousef lag , xymouse , $ 


; Image 1 and its dimensions. 

; Image 2 and its dimensions. 

;Coordinates of draw windows (LL corner) . 
;Current magnification of images 1 and 2. 
;Windov IDs for display widgets. 

; Window IDs for detail displays. 

;Size of the draw windows. 

;Qutput file flag and unit number. 

;Flag and position of mouse. 
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xl, yl, xlflag, 
x2, y2, x2flag, 
xr, xpr , 
detsz 

COMMON vidget.ids, 
irbase , 
drawl, draw2, 
detaill, detail2, 
textlx, textly, 
text2x, text2y 

; Constants : 
d2 = detsz/4 

xpad = 0 
ypad = 0 

CASE vid of 


$ ;Selected points in image 1 and flag. 
$ ;ditto for image 2. 

$ ; Returned arrays. 

;Size of detail windows 

$ 

$ ;Main base ID. 

$ ;Drawing areas. 

$ ;Detail view areas. 

$ ;Text output areas for image 1. 

;ditto for image 2. 


detaill: begin 

winid = detlid 
d2 = d2/magnl 
d = 2*d2 

subim = BYTARR(d,d) 
xlo = (xm-d2+l ) 
xhi = (xm+d2) 

if xlo It 0 then xpad = d-xhi-1 
ylo = (ym-d2+l) 
yhi = (ym+d2) 

if ylo It 0 then ypad = d-yhi-1 

subim(xpad , ypad) = BYTSCL( iml [(xlo>0) : (xhi<(iml_xO-l) ) . (ylo>0) : (yhi< ( lal yO-1))] ) 

end 


detail2: begin 

winid = det2id 
d2 = d2/magn2 
d = 2*d2 


subim = BYTARR(d.d) 
xlo = (xm-d2+l) 
xhi = (xm+d2) 


if xlo It 0 then xpad = d-xhi-1 
ylo = (ym-d2+l) 
yhi = (ym+d2) 

if ylo It 0 then ypad = d-yhi-1 

subim(xpad,ypad) = BYTSCL( im2[(xlo>0) : (xhi<(im2_x0-l)) , 


(ylo>0) : (yhi < (im2_y0-l ) )] 


) 
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end 

end 

VSET,vinid 

TVSCL , CONGRID (subim , detsz , detsz) 

TVCIRCLE , 20 ( detsz/2 , detsz/2 , color=200 

PLOTS, [0,detsz-l] , [detsz/2, detsz/2] ,color=200, /device 

PLOTS, [detsz/2, detsz/2] , [0,detsz-l] .color =200, /device 

RETURN 

END 


PRO Image.Register, imagel, image2, zret, xpret, $ 

CENTER=xycen, OUTPUT=output 

, MAIN PROGRAM: Creates the main vindov widget and registers it with XManager. 


COMMON im_register, $ 

ial, iml_xO, iml.yO, $ 

ia2, im2_x0, im2_yG, $ 

xyvisl , xyvis2, $ 

nagnl, magn2, $ 

drawl id ,draw2id , % 

detlid, det2id, $ 

winsz, $ 

out file ,unitnum, $ 

nous ef lag ,xymouse , $ 

xl, yl, xlflag, $ 

x2 , y2, x2f lag , $ 

xr, xpr , $ 

detsz 

COMMON widget. ids, $ 

irbase, $ 

drawl, draw2 , $ 

detain, detail2, $ 

text lx , text ly , $ 

text2x , text2y 


;0N_ERROR f 2 


; Image 1 and its dimensions. 

; Image 2 and its dimensions. 

;Coordinates of draw windows (LL corner). 
;Current magnification of images 1 and 2. 
;Vindow IDs for display widgets. 

; Window IDs for detail displays. 

;Size of the draw windows. 

;0utput file flag and unit number. 

.* Flag and position of mouse. 

; Selected points in image 1 and flag. 
;ditto for image 2. 

; Returned arrays. 

; Size of detail windows 


;Main base ID. 

;Drawing areas. 

;Detail view areas. 

I Text output areas for image 1. 
;ditto for image 2. 


check for previous registration: because of 


common block, can only have one 
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;realization of XYVieu operating at a time. 

if (Xregistered('Image.Register') ne 0) then MESSAGE, 'Widget is already registered' 

if N.ELEMENTS ( image 1 ) eq 0 then iml=DIST(512,512) else iml = imagel 
if N.ELEMENTS (image2) eq 0 then im2=DIST(512,512) else in>2 = imaee2 
szl = SIZE(iml) 6 

sz2 = SIZE(im2) 

if szl[0] ne 2 or sz2[0] ne 2 then MESSAGE, 'One or both of the input images is not 2 

; Initialize output state: 
outfile = 0 

if KEYWORD.SET (output) then begin 
outfile = 1 

OPENV , unitniun , STRING (output) , /GET.LUN 
end 


; Initialize the selected points array (l-D coordinates)- 
xl = [-1.) 
x2 = xl 

yi = [-U 
y2 = yl 
xlflag=0 
x2flag=l 

mouseflag=0 


;Widget creation: 

> 1 • Main base: 

lrbase - WIDGET_BASE( TITLE=’Iaage Registration'. MBAR=bar ) 
pdstruct = {flags:0,name: "} 
desc = REPLICATE (pdstruct, 17) 

desc . f lags = [l ,0,0,2 , 1,2. 1 , 1 .0.0.2, 1 .0 2 2 1 21 

aesc.name=[ 'File 1 , $ 

'Save', 'Save As', 'Exit* % 

'Edit 1 , $ 

'Undo Selection' ,$ 

'Display', $ 


'Selections', $ 

'Print', 'Circle', 'Highlight 

'Color Table' , $ 

'Xloadct . . . ' , 'XPalette . . . ’ " $ 

'Transform’,? 

’Affine’] 


aainmenu = CV_PDME.VU( bar, desc. 


/RETURN. FULL.NAME, UVALUE= 


’menu’, /MBAR ) 


: returning 
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;2. Display areas: main windows are 480x480 pixels, detail windows are 256x256 
minsz = 380 


winsz = 480 
detsz = 256 


iml-XO = szl[l] 
iml.yO = szl [2] 
magnl = 1 

imls = MIN([iml_xO,iml_yO] .xymini) 


im2_r0 = sz2 [1] 
im2_y0 = sz2[2] 
magn2 = 1 

im2s = MIN ( [im2_x0 , im2_y0] , xymin2) 
iain = MIN(imls, ia2s) 

if imin It minsz then winsz = minsz else minsz = winsz 


drawbase = WIDGET_BASE( irbase, GROUP_LEADER=irbase, /FRAME, /ROW ) 

desc2 = REPLICATE(pdstmct , 11 ) 
desc2. flags = [1 ,0,0,0, 0,0, 2, 1,0, 0,2] 
desc2.name = ["Zoom",$ 

# 1:4V1:2V1:1 V2:l V4:1V8:1" # $ 

"Move ' , $ 

"Image Center" , "House. . . * , "Enter. . . *J 
; Image 1 windows and controls 

/col™, 

] WIDCET. LABEL ( drawl.base, VALUE* 'Reference Image' ) 

drawl = WIDGET_DRAW( drawl.base, RETAIN=2, CRAPHICS.LEVELM . $ 

XSIZE=iml_xO, YSIZE=iml_yO, FRAME=2,$ 

/BUTTON. EVENTS , /MOTION.EVENTS , UVALUE = 'drawl'. $ 

lowl = WIDGET BASEf / S Tk L ' * 0LL - SI ZE=winsz , Y_SCROLL_SIZE=winsz) 

con, , WIDCEt'bacfi ' / ' UC “- LEFT - ««0»PJ.EUia.irta,.,/IU)w) 

WIDGET.BASEC lowl, CROUP. LEADER=irbase /COL) 

-Ux^^DG™" 1 * ? e3C2> /R£TO « L -^. ™^'men U , ) 
tlx WIDCET_TEXT( coni, FRAME=0, SCR_XSIZE=100 , UVALUE= • text lx • ) 

extly = WIDGET_TEXT( coni, FRAMED, SCR.XSIZE=100, UVALUE= ' text l v ■ 

e-aill - WIDGET. DRAW ( lowl, RETAIN=2 , CRAPHICS.LEV=1 , $ ' 

SCR_XSIZE=256 , SCR_YSIZE=256 , FRAME=2, $ 

/BUTTON. EVENTS, UVALUE= 'detl ' ) 

; image 2 windows and controls 

r^^DGET UBEuT E( draaba3e ’ G * 0UP -^DER=drawbase. /COLUMN) 
jun* WIDGET.LABEL( draw 2 _base, VALUE* ’Test Image’ ) 

WIDGET_DRAW( draw2.base, RETAIN=2, CRAPHICS LEVEL =1 J 

XSIZE=im2_x0, YSIZE=im2_yO, FRAME=2, $’ 


NASW-98008 


NASW-98008 


Final Report: Spring 2001 


36 


4.3. IDL Destretching Tools 


FUNCTION destretch, imO, iml, gridnest, fvfac, grid, stretch_clip,$ 
OVERLAP =over lap 


;Uses ANA routines GRIDMATCH and STRETCH (available via DLMs 
; at http://ana.lmsal.com/anaidl) to destretch IH1 onto IMO. Useful only 
; for destretching a single pair of images. Time series require 
; special grid handling and smoothing. 

;GRIDNEST is an N-dimensional integer array in vhich 
; GRIDNEST(n) = the number of subtiles to be used on the 
; nth destretching iteration. The tiles have the same aspect 
; ratio as the images. Thus, if the image is not square, the 
; user must ensure that the grid size is sufficient in longest 
; dimension. 

; FVFAC: apodization gaussian FVHM = least tile width • fvfac 

; STRETCH _CLI P : array: mar i m u m displacement allowed, in pixels 
; for each grid in gridnest. 

;0UTPUT: IMIS = destretched IM1 

* GRID = FLTARR (2, MAX (gridnest ), MAX (gridnest) ) 

» Note that this grid is only the final incremental 

J grid applied to the image, not the full dis- 

» placement grid needed to go from iml to imO! 

0N_ERRQR , 2 

szO = SIZE(imO) 
szl = SIZE(iml) 

if szO(O) ne 2 or szl(O) ne 2 then begin 
MESSAGE, 1 Images must be 2-D' 

RETURN 

end 

if szO(l) ne szl ( 1 ) or sz0(2) ne szl(2) then begin 
MESSAGE, ’Images must be same size’ 

RETURN 

end 
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if N.ELEMENTS(stretch.clip) eq 0 then str=REPLICATE (20 , N_ELEMENTS (gridnes t ) ) else $ 
str=stretch_clip 

;default tile overlapping = 10X 

if KEYWORD. SET (overlap) then ovr = overlap else ovr = 0.10 
imls = iml 

for i=0,N_ELEMENTS(gridnest)-l do begin 

ni = gridnest(i) 
ngx = sz0(l)/ni 
ngy = sz0(2)/ni 

grid.x = (IKDGEN(ni)*ngx ♦ ngx/2)# (INTARR(ni)+l) 
grid.y = (INTARR(ni)+l)#(IHDCEN(ni) *ngy + ngy/2) 

vind_x = FIX (ngx+ovr*ngx) 

if wind.x mod 2 eq 0 then yind_x = vind_x ♦ 1 
wind.y = FIX(ngy+ovr*ngy) 

if wind.y mod 2 eq 0 then vind.y = wind y 4 1 

fvhm = FIX(MIN([ngx*fyfac,ngy*fyfac])) 

grid = GRIDHATCH(imO,ials,grid_x,grid_y ,vind_x,wind_y,fvhm,Btr (i)) 

imls = STRETCH(imls ,grid) 


end 

RETITRN, imls 
E?fD 


PRO ds.series 


lntmpl , i0, il, outmpl , grids, clips, sfac, delta, $ 
MAXDISP=Bardisp, LOGFILE=logf ile , SILEMT=silent 


szg = SIZE(grids) 
ngrid = szg(l) 

°arg = grids (ngrid-1) 

.Default: no unsharp masking 
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if N_ELEMENTS(sfac) eq 0 then sfac = il-iO 

;The maximum allowed displacement value in pixels: 
if KEYWORD.SET (maxdisp) then m&xstr = maxdisp else maxstr=20 
;Shutup if asked 

if KEYWORD. SET (si lent) then loud = 0 else loud = 1 
; Redirect output to a logfile if asked 

;if KEYWORD.SET (logfile) then SPAWN, 1 ... .redirect stdout to file...* 


tO = SYSTIME(l) 

if loud then PRINT, Calculating displacement grid...' 

FZREAD,ref im,fns(intmpl, iO) ,h ;Note the use of the RIG2 files (rigid aligned) 

FCWRITE,ref im,fns(outmpl,iO) ,h,5 
delta « FLTARR(il-iO+l ,2,maxg,maxg) 
for i=iO+l,il do begin 

tOi = SYSTIHE(l) 

if loud then PRINT, ’Working on image ’ ,i 
FZREAD, im,fns(intmpl,i) ,h 
dq = dsgridnest (ref im,im,gTids, clips) 
badg = WHERE(dq gt maxstr, num) 
if num gt 0 then begin 

if loud then PRINT, 1 Accumulated offsets > maxdisp: *,mim 
dq(badg) =0.0 

end 

delta(i-iO, *,*,•) = dq 

if loud then PRINT, ’Mean dx, dy = ’ ,STRTRIM(AVC(ABS(dq(0, • , •) ) ) 2) ’ ' $ 
STRTRIM(AVG(ABS(dq(l, •,•))) , 2 ) 

if loud then PRINT, 'Max displ = ’ ,STRTRIM(MAX(ABS(dq)) ,2) 


refim = im 
im = 0 

if loud then begin 

PRINT, 'Time for image ’ ,STRTRIH(i ,2) , ’ = ’ ,STRTRIM(SYSTIHE(l)-tOi ,2) , ' seconds’ 
PRINT,’ ’ 

end 


end 

if loud then PRINT, 'Time for displacement grid calculation: ’ ,SYSTIM£( 1 ) -tO , ' seconds 
;Zero out the IEEE NaN values from GRIDMATCH: 
delta (WHERE (FINITE (delta) eq 0)) = 0.0 

; Linear detrend and unsharp the displacements: 
delx = REFORM (delta (* ,0, * , *) ) 
dely = REFORM(delta(* , 1 , * , *) ) 
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dsz = SIZE(delx) 
nt = dsz(l) 
xvec = INDGEN(nt) 
ngx = dsz (2) 
ngy = dsz (3) 
delta = 0 

if loud then PRINT, ’Detrending and unsharp masking the displacements...' 
tO = SYSTIME(l) 

for j=0,ngy-l do for i=0,ngx-l do begin 

xq = TOTAL(delx(*,i,j) ,/cuau) 
cf = POLY_FIT(xvec,xq, 1 ,yf it) 
xq = xq - yfit 

xq = xq - SMOOTH (xq.sfac, /EDGE) 

delx(0,i, j) = xq - xq(0) ; insure good starting point. 

xq = TOTAL(dely(*,i,j),/cuau) 
cf = POLY_FIT(xvec,xq,l,yfit) 
xq = xq - yfit 

xq = xq - SMOOTH (xq, sfac , /EDGE) 
dely(0,i,j) = xq - xq(0) 


end 

if loud then PRINT. 'Detrend and smooth time = ’ ,SYSTIME(l)-tO. ' seconds’ 

delta = FLTARR(nt ,2, ngx, ngy) 
delta(*,0,*,*) = delx 
delta(* , 1 , *,*) = dely 
delx = 0 
dely = 0 

;Apply the grid: 
tO = SYSTIME(l) 

if loud then PRINT, Applying processed grids to images... 1 
for i=iO+l , il do begin 

FZREAD, im, fns(intmpl ,i) ,h 
gridi = REFORM (delta (i-iO ,*,*,»)) 
gstd = STDEV (gridi , gavg) 

bads = VfHERECABS (gridi) gt gavg + 3»g 3 td,nbad) 
if nbad gt 0 then gridi(bads) = 0 
im = STRETCH (im , gridi ) 

fileout = fns(outmpl,i) 

if loud then PRINT , 'Writing \fileout 
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FCWRITE, im,f ileout ,h,5 


if loud then PRINT, 'Time for grid application: ' ,SYSTIME(l)-tO, * seconds' 
if loud then begin 
PRINT , * » 

PRINT,'Done destretching the tine series. 1 


RETURN 

END 


4.4. IDL FFT and P-mode Filtration Routines 


•include <stdio.h> 

•include <math.h> 

•include "export. h" 

•include "nr.h” 

•include "nrutil.h" 

•define ARCS int, IDL.VPTRf] 

•define ARRLEN(arr) (sizeof (arr)/sizeof (arr[0])) 

•define SVAP(a.b) f teap*(a) ; (a)«(b) ; (b)=f temp 
•define NRANSI 

static void IDL_fourfs( int argc, IDL.VPTR »argv) ; 

void fourev(FILE .file[5], int ma. int -nb, int int .nd) 

void print_vptr( IDL.VPTR, char « ); 

int IDL.Load (void) 

static IDL.SYSFUN.DEF2 nev.pro[] = { 

{ (IDL.FUN.RET) IDL.fourfs, "FOURFS”, 3, 3, 0, 0> 


return IDL_SysRtnAdd(nev_pro, FALSE, ARRLEN(nev.pro) ) ; 



static void IDL.fourfs ( int argc, IDL.VPTR .argv) 


/ 


Singleton’s algorithm for FFT 


of a dataset stored on disk. 
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Modified slightly from Numerical Recipes in C, 2nd Edition, 
Section 12.6. 

NOTE! ! Data must be COMPLEX valued. For each read, the 
routine assumes it’s reading a (real,imag) value pair 
stored in native floating point in the disk files. Thus 
for a (nx=128,ny=128,nt=256) image cube, the resulting 
disk file A containing 1/2 of the data would be 
128L*128*256*sizeof (float) *2 bytes in size. That extra 
little factor of 2 caused me no end of pain when first 
trying to use this routine on real valued data. 


Inputs: 

fnames [] 

string array, contains the filenames 
of the 4 storage arrays on disk. 


nn[] 

int array, contains the length of each 
dimension of the data array. 


isign 

int, forward or reverse transform flag 
1 = forward transform 
•1 = inverse transform 

Outputs: 

none. 



chan: warn; 

int cc,na,nb,nc,nd; 

float tempr , tempi , *af a, *afb, *af c 

double wr,vi,vpr,vpi,vtemp, theta 

static int mate[5] = {0,2, 1 , 4 , 3} 

FILE *f ile [5] ; 


IDL.FILE.STAT fstat; 

IDLJJLONG i , j » jl2, jk,k,kk f n=l t mm,kc=0,kd,ks,kr 
IDL.VPTR fnames ,nnv,dum[2] ; 

IDL_L0NG luns [4] , *nn ; 

IDL_VARIABLE unit, name; 

IDL.STRING *fn; 


nr ,ns ,nv; 


int ndim, isign; 


^define KBF 1048576 

/•Setting KBF to 128 or 512 gave totally screwed up transform 

on the test case run on my SGI 02. KBF=1024 worked, but no faster than 256 

Leaving it at 256... •/ 

gger test cases take 25 hours to transform! Need to set KBF MUCH large 


?er 
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in order to make the routine reasonable - or look into virtual memory tricks...*/ 

/* printf ("argc = Xd\n",argc) ; •/ 

/•get file pointers set up based on argv[0] array*/ 
fnames= argv[0]; 

/*print_vptr(luns,"luns"); •/ 

IDL.ENSURE.STRING (f names) ; 

IDL.ENSURE. ARRAY (f names) ; 
if (fnames->value.arr->n_dim != 1) 

IDL.Message (IDL.M.NAHED.GENERIC , IDL.MSC.LONGJMP , 

"Logical unit argument must be a 1-D array."); 
if (fname3->value.arr->n_elts != 4) 

IDL.Message (IDL.M.NAHED.GENERIC , IDL.MSC.LONGJHP , 

"Wrong number of filenames passed."); 

/• printf("\nf names type: ’/.d\n",fnames->type) ; •/ 

/• printf ("fnames nelt = Xd\n",fnames->value.arr->n_elts) ; •/ 

/•get filenames */ 
printf ("\nStorage Files :\n"); 
fn = ftfnames->value,arr->data[ 0 ]; 
for(i=0;i<4;i++) •( 

/• printf ("filename Xd length » */.d\n\i ,fn->slen) ; •/ 

/• printf ("filename Xd type * Xd\n",i,fn->«type) ; •/ 

printf ("File Xd = Xs\n",i+l,fn->s) ; 
fn++; 

} 

/• Open files •/ 

printf ("\nOpening f iles . . . \n") ; 

fn = &fnames->value.arr->data[ 0 ] ; 

f or (i=0 ; i<4 ; i++) { 

/* get file unit */ 
dum[0] = tunit; 
unit. type = IDL.TYP.LONG; 
unit. flags = 0; 

IDL.FileGetUnit ( 1 ,dum) ; 
luns[i] = dum[0] ->value . 1 ; 

/* set up filename •/ 

/* printf ("Opening LUN Xd\n".(int) luns[i)) ; ./ 
name. type = IDL.TYP.STRING ; 
najne. flags = IDL_V_CONST; 
name .value .str. s = fn->s ; 
name. value. str.slen = sizecf (fn->s) -1 ; 
name . value . st r . stype = 0; 
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dum[l] = fcname; 

/* NOTE: Extern. Dev. Guide function prototype incorrect: needs one more arg */ 
IDL_File0pen(2, dum, (char *) 0, IDL.OPEN.R | IDL.OPEN.W, 

IDL_F_UNIX_F77 , 1,0); 

fn++; 

> 

printf ("Files successfully opened. \n"); 

/♦Check file status */ 

printf ("\nEnsuring file status. . An") ; 

for(i=0;i<4;i++){ 

if (IDL_FileEnsureStatus(IDL_HSG_LONGJMP, (int) luns[i], 

IDL.EFS.OPEN | IDL.EFS.READ I IDL.EFS.VRITE) ) 
printfC'File '/.d is open for reading and writing. \n”, 

(int) luns[i]); 

> 

printf ("Status ensured. \n") ; 

printf ("\nGetting FILE pointers. . An”) ; 
for (i=0;i<4;){ 

IDL_FileStat((int) lunsti], fcf stat) ; 

/• printf ("Filename from FileStatO : */,s\n" ,f stat .name) ; •/ 

f ile[++i] = fstat.fptr; /»Hote unity offset required by NR! ! •/ 

> 

printf ("FileStats complete. \n") ; 

/» Establish dimensions of data array ♦ / 
nnv = argv[l] ; 

/•print_vptr(nnv,"nnv") ; ♦ / 

IDL_ENSURE_SIMPLE(rmv) ; 

IDL.ENSURE. ARRAY (nnv) ; 

nn = (IDL.LONG ») nnv-> value. arr->data; 

if (nnv-> value. arr->n.dim != 1 ) 

IDL_Message(IDL_M_NAMED_CEXERIC , IDL_HSG_L0NG JMP , 

"Array dimensions argument mu 3 t be a 1-D 2 irray."); 
ndim = nnv->value . arr->n_elts ; /♦confusing name, but keep NR names intact... ♦/ 

/» get isign •/ 

isign = argv [2] -> value . i ; 

/» printf ("isign= , /.d\n", isign) ; ♦/ 

if (isign == 1) { 

for (i=0;i<ndim;i++) { 
if (nn[i] < KBF) warn=l; 

> 

if (warn) IDL_Message(IDL_K_NAMED_GENERIC , IDL.HSG.INFO , 
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"At least one array dimension < KBF record length: Transform will be wrappi 

> 

nn -= 1; /* unity offset for all NR routines!! */ 

/* for (i=l ; i<=ndim; i++) printf ("nn('/,d) = '/.d ",i,nn[i]); •/ 

/* Here begins the NR routine (almost) unmodified from the original: •/ 

/• printf ("\n\nsizeof (float) = %d\n",sizeof (float)) ; •/ 
afa=vector(l ,KBF) ; 
afb=vector (1 ,KBF) ; 
af c=vector(l ,KBF) ; 
for (j=l; j<=ndim; j++) { 
n *= nn[j] ; 

if (nn [j] <= 1) IDL.HessagedDL.H.NAMED.CENERIC, IDL.MSC.LONGJMP , 

"Invalid float or wrong ndim in fourfs"); 

> 

/• printf ("\nn = '/.d\n' , ,n); •/ 
nv=l ; 

jk=nn [nv] ; 
nm-n ; 
ns=n/KBF; 

/• printf ( M ns * %d\n N ,ns); •/ 
nr=ns >> 1; 

/* printf ("nr = '/.dXn'Snr) ; •/ 

kd=KBF » 1; 
ks=n; 

fourev(f ile,fcna, fcnb,&nc,fcnd) ; 

/* First phase of the transform starts here •/ 

/• printf ("\nStarting first phase ... \n") ; •/ 

/• printf ("mm = '/.d\n" ,mm) ; •/ 

/• cnt=0; */ 

for (;;) { /*Start of the computing pass. •/ 

/• Done M-K-l times, where n = 2‘M and KBF=2‘K •/ 

/• printf ("count = 7.d\n" ,++cnt) ; •/ 
theta=isign*3. 141592653589793/ (n/mm) ; 
vtemp=sin(0.5*theta) ; 
wpr = -2 .0*wtemp*wtemp; 
vpi=sin(theta) ; 
wr=l . 0 ; 
vi=0 . 0 ; 
mm >>= 1 ; 

/• printf (" mm = '/.d\n" ,=a) ; »/ 
for ( jl2=l ; j 12<=2; jl2++) { 

/* printf (” j 12 = 7.d\n" , j 12) ; */ 

kr=0 ; 
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/* printf ("ftell(file[na]) : %d\n" ,f tell(f ile[na] )) ; ♦/ 

/* printf ("f tell (file [nb]) : 7*d\n” ,f tell (file [nb] )) ; •/ 


b 


do { 


printf ("kr = 7.d\n M ,kr); */ 

cc=fread(Jtaf a[l] , sizeof (float) ,KBF,f ile [na] ) ; 
if (cc != KBF) { 

printf ("KBF = Xd, cc = Xd\n" f KBF f cc) ; 
IDL.Hessage(IDL.M.NAMED„GENERIC f IDL_MSG_L0NG JMP , 

"read error in fourfs 1"); 


> 


cc=f read(fcafb[l] t sizeof (float) ,KBF,f ile [nb] ) ; 
if (cc ! = KBF) { 

printf ("KBF = Xd, cc = XdWjCBF.cc) ; 

IDL.Message (IDL.H.NAMED.GENERIC , IDL.MSG.LONC JHP , 

"read error in fourfs 2"); 


} 


for (j=l ; j<=KBF; j+=2) { 

tempr=((f loat)irr)*aIb[j] - ((float )vi)*afb[j + l] ; 

tempi=((f loat)ui)*afb[j] ♦( (float )vr)*afb[j + l] ; 

afb[j]=afa[j]-tempr ; 

afa[j) += tempr; 

afb[j + l]=af a[j + l]-tempi; 

afa[j+l] += tempi; 

} 

kc +« kd; 
if (kc == mm) { 
kc=0; 

vr=(vtemp=vr) *vpr-vi«vpi+vr ; 
vi=wi*vpr+vtemp*vpi+vi ; 


cc=fwrite(*afa[l] .sizeof (float) , KBF, file[nc]) ; 
if (cc ! = KBF) { 

printf ("KBF = 7.d, cc = 7,d\n" .KBF.cc) ; 
IDL_He3sage(IDL.M_NAMED. GENERIC , IDL.HSG.LONGJMP , 
"vrite error in fourfs 1"); 


> 


cc=fwrite(4afb[l] .sizeof (float) , KBF, file[nd]) ; 

if (cc ! = KBF) { 

printf ("KBF = 7,d, cc = */.d\n" .KBF.cc) ; 
IDL_Message(IDL_H_NAHED_GENERIC, IDL.HSG.LONGJMP, 

"write error in fourfs 2"); 
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} while (++kr < nr); /*end of the big DO loop •/ 

if (j 12 == 1 44 ks != n U ks == KBF) { 
na=mate[na] ; 
nb=na; 

> 

if (nr == 0) break; 

> 

f ourew(f ile ,4na,4nb ,4nc ,4nd) ; /*Start of the permutation pass • / 

jk »= 1; 
while (jk == 1) { 
mm=n; 

jk=nn[++nv] ; 

> 

ks >>= 1; 
if (ks > KBF) { 
for (jl2=l; jl2<=2; jl2++) { 

for (kr=l ;kr<=ns;kr+=ks/KBF) { 
for (k=l ; k<=ks ;k+=KBF) { 
cc=fread(4afa[l] .sizeof (float) , KBF, file[na] ) ; 
if (cc ! = KBF) IDL.Hessage (IDL.M.NAHED.CENERIC , IDL.HSC.LONCJMP, 

"read error in fourfa 3"); 
cc=fwrite(iaf a[l] .sizeof (float) ,KBF,f ilefnc] ) ; 
if (cc ! = KBF) IDL.Hessage ( IDL_H_NAHED_GENERIC , IDL.HSC.LONCJMP, 

"read error in fourfa 4"); 

nc=mate[nc] ; 

> 

na=mate [na] ; 

> 

fourew(f ile,4na,4nb,knc,bid) ; 

> else if (ks == KBF) nb=na; 
else break; 

> 

pr intf ( "Start ing second phase ... \n") ; •/ 

/•Second phase starts here. Done K times where KBF = 2‘K •/ 

j=i; 

for (; ;) { 

theta=isign*3 . 141592653589793/ (n/mm) ; 

wtemp=sin(0 . 5* theta) ; 

vpr = -2 .0*vtemp*wtemp; 

vpi=sin (theta) ; 

wr=l .0 ; 

vi=0 . 0 ; 
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mm >>= i ; 

/* printf(" mm = ‘/d\n" ,mm) ; */ 
ks=kd; 
kd »= 1; 

for (jl2=l; jl2<=2; jl2++) { 
for (kr=l;kr<=ns;kr++) { 

cc=f read(4af c [1] .sizeof (float) ,KBF,f ile[na] ) ; 

if (cc != KBF) IDL.Message (IDL_M_NAMED_GENERIC , IDL.MSG.LONGJHP, 

"read error in fourfs 5"); 

kk=l; 
k=ks+l ; 
for (;;) { 

tempr=((f loat)wr)»afc[kk+ks] -( (float )vi)*afc [kk+ks+l] ; 
tempi=( (float) vi) *afc[kk+ks] + ( (float )wr)»afc[kk*ks+l] ; 
afa[j]=afc[kk]+tempr; 
afb[j]=afc[kk]-tempr; 
af a[++j]=af c [++kk)+ tempi ; 
afb(j++] =af c [kk++] -tempi ; 
if (kk < k) continue; 
kc 4= kd; 
if (kc == mm) { 
kc=0; 

wr = ( ut emp =wr ) • wpr- wi • up i 4wr ; 
wi=wi*wpr+wtemp»Bpi4Hi ; 

> 

kk 4= ks; 

if (kk > KBF) break; 
else k=kk+ks ; 

> 

if (j > KBF) { 

cc=f write (4af a[l] .sizeof (float) ,KBF,f ilefnc] ) ; 

if (cc ! = KBF) IDL.Kessage ( IDL_H_NAHED_GENERIC , 1DL.MSG.L0NCJMP , 

"write error in fourfs 3 M ); 
cc=fwrite(4af b [1] .sizeof (float) ,KBF,f ile[nd] ) ; 

if (cc ! = KBF) IDL.Hessage ( IDL_M_NAMED_GENERIC , IDL.KSG.LONCJMP , 

"write error in fourfs 4"); 

> 

} /* end of kr loop */ 
na=raate [na] ; 

} /* end of j 12 loop */ 

f ourew (f ile , fcna , &nb , 4nc , tnd) ; 
jk »= 1; 

• printf ("jk = 7.d\n', jk) ; ./ 
if (jk > 1) continue; 
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mm=n; 
do { 

if (nv < ndim) { 
jk=nn[++nv] ; 

/* printf ("remainders: jk = Xd\n M ,jk); ♦/ 

} 

else { 

free.vector (af c, 1 ,KBF) ; 
free_vector(afb,l,KBF) ; 
free_vector(afa f l,KBF) ; 

/♦Free the file units*/ 

printf ("\nFreeing the file units ... \n") ; 

for (i=0;i<4; i++) { 

/• printf ("Freeing LUN Xd\n", (int) luns[i]); •/ 

/• setup IDL.VPTR •/ 
unit. type = IDL.TYP.LONC ; 
unit. flags = 0; 
unit. value. 1 = luns[i]; 
dum[0] = fcunit; 

IDL.FileFreeUnit (1 , dua) ; 

> 

return; 

> 

> while (jk == 1); 

} 

> 

Sundef KBF 
tfundef NRANSI 


void fourew(FILE •file[5], int .na, int . n b, ant int -nd) 


int i ; 

FILE *ftemp; 

for (i=l ; i<=4 ; i++) revind(f ile [i] ) ; 
SWAP (f ile [2] , file [4] ) ; 

SWAP (file [l] .file [3] ) ; 

*na=3 ; 

*nb=4 ; 

*nc=l ; 

♦rid=2 ; 
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•undef SWAP 


/* Applies sub-v or sub-fundamental space-time filters in 
Fourier domain using the temporary file structure given by 
fourfs.c. Useful for p-mode filtering large datasets that 
don't fit in RAM. */ 

•include <stdio.h> 

•include <math.h> 

•include "export. h" 

•include "nr.h" 

•include "nrutil.h" 

•define ARRLEN(arr) (sizeof (arr)/sizeof (arr[0])) 

•define SWAP(a.b) f temp=(a) ; (a)=(b) ; (b)=f temp 
•define NRANSI 


IDL.VPTR IDL.pmodef ilter ( int argc, IDL.VPTR «argv ); 
void filesetup( IDL.VPTR ); 

int IDL.Load(void) 

{ 

static IDL.SYSFUN.DEF2 nev.pro[] = { 

{ (IDL.FUN.RET) IDL.pmodef ilter , "PMODEFILTER” , 5, 5, 0, 0) 


return IDL_SysRtnAdd(nev_pro, TRUE, ARRLEH (nev.pro) ) ; 


•define NF 2 
FILE *f ile [NF] ; 
IDL.LONG luns [NF] 


/•current algorithm stores transformed data 


in 2 storage arrays •/ 




IDL.VPTR IDL.pmodef ilter ( int argc, IDL.VPTR *argv ) 


Inputs: f names [2] string array: filenames of temporary files. 

[nx.ny.nt] : dimensions of ORIGINAL (ie. not transposed) data cube, 
dxy: spatial scale units in km/pixel. 
dt . temporal scale units in seconds, 
vtop: cutoff velocity. 
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Outputs: 3/7/01. last image read from data cube. 

see pmf _test.pro for testing routine. 

****************************************** ••••*••*•*•*••**••• 

char *dataO,*datal ; 

long *xyt,i, j ,k, jk, jkc.ll; 

unsigned int ndim.f sk,sf=sizeof (float) ; 

float dk , dkx , dky , dt ,dv,dxy,kx,ky,kc,v,vtop,w t vc; 

size_t cc; 

IDL. VARIABLE unit; 

IDL.VPTR fn,nnv,dum[2] ,tmpO,tmpl ,dxyv,dtv, vtopv; 

IDL_L0NG nc,nf ,nt,nt2,nx,nx2,ny,ny2,ntny[2] ; 


/•open files and assign logical units •/ 
f n = argv [0] ; 
f ilesetupC fn ) ; 


/•Get dimensions of original data cube •/ 
nnv = argv[l] ; 

IDL_ENSURE_SIMPLE(nnv) ; 

IDL_ENSURE_ ARRAY (nnv) ; 
if (nnv->value . arr->n_elts != 3) 

IDL.Message ( IDL.M.NAMED.CEJfERIC , IDL.MSG.LONC JMP , 

"Data array must be 3-D."); 
xyt = nnv->value . arr->data; 

nx = xyt [0] ; 
nx2 = nx/2; 
uy = xyt [1] ; 
ny2 = ny/2 ; 
nt = xyt [2] ; 
nt2 = nt/2 ; 

nc = nfny; /• number of complex elements in each u-ky image./ 

nf 2 * nC ; number float elements in each v-ky image (real, mag) •/ 

printf ("Dimension of images in series from nn array: \n"); 
printf("\tnx = ’/,d\n",nx); 
printf ("\tny = 7.d\n",ny); 

printf ("Number of images in series from nn array: \n">- 
printf ("\tnt = 7.d\n M ,nt); 

printf ("Number of floating point variables in a w-ky image- \ n ") • 
printf ("\tnf = */.d\ n ",nf); 8 ' 
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/•Units and cutoff velocity - note assumptions on data types!*/ 
dxyv = argv [2] ; 
if (dxyv->type != 4) 

IDL.Message (IDL_M_NAHED_GENERIC , IDL.MSG.LONGJMP, 

"DXY argument must be floating point."); 
dry = dxyv->value . f ; 

dtv = argv [3) ; 
if (dtv->type != 4) 

IDL.Message (IDL.M.NAKED.GENERIC, IDL.MSG.LONGJMP, 

"DT argument must be floating point."); 
dt = dtv->value.f ; 


vtopv = argv [4] ; 
if (vtopv->type != 4) 

IDL_Message(IDL_M_NAMED_CEVERIC, IDL.HSG.LONCJHP, 

"VTOP argument must be floating point."); 
vtop = vtopv->value.f ; 


/•Fourier domain units •/ 

dv = l./nt/dt; 

vc = l./2/dt; 

dkx = l./nx/dxy; 

dky = l./ny/dry; 

kc = l./2/dry; 

printf ("dv = 7f\n”,dv); 

printf ("dkx = ‘/.f \n" ,dkx) ; 

printf ("dky = 7.f\n",dky); 


/• Get temporary array for v-k 7 image planes •/ 
ntny [0] =nt ; 
ntny [1] =ny; 

dataO = IDL.MakeTempArray (IDL.TYP.COMPLEX , 2 , ntny , IDL. ARR. INI _N0P , ttmpO) ; 
datal = IDL.MakeTempArray ( IDL.TYP.COMPLEX , 2 , ntny , IDL_ ARR_ IN I _N0P , ttmp 1 ) ; 

/• Loop through storage arrays. Dimensions of each storage file = ( n t ny nx/2) 
Note symmetry: filel is mirror of fileO in time axis for real input*/’ 

/• 3/7/01: This is now successfully reading the v-ky TRANSPOSE imago 
from the data cube. NEXT: figure out how to read, 

filter it, and transpose it back so it ’ 3 ready for the IFFT •/ 
for (i=0;i<nx2;i++) { 

/♦read NxM complex iinage •/ 

/• read a v-ky slice image from 1st storage file */ 
fsk = fseek(file[0] ,i»nf»sizeof (float) .SEEK. SET) ; 

cc = f read(4data0[0] .sizeof (float) ,nf ,file[ 0 ] ) ; 
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if (cc ! = nf) IDL.MessagedDL.M.NAHED.GENERIC, IDL.MSG.LONGJMP , 

"Read failure on file[0]."); 

/♦ read a w-ky slice image from 2nd storage file which is in reverse kx order*/ 
fsk = f seek(file[l] , (nx2-i-l)*nf*sizeof (float) ,SEEK_SET) ; 
cc = fread(kdatal[0] ,sizeof (float) ,nf,file[l]) ; 
if (cc != nf) IDL.Mes sage (IDL.M.NAHED.GENERIC t IDL.MSG.LONGJHP, 

"Read failure on f ile [1] . ,f ) ; 

/♦condition the transform images */ 
kx = i*dkx; 

for ( jk=0; jk<=(sf *nf ) ; jk++){ /*loop over all bytes in the image array •/ 

jkc = jk/(2*sf); /*pixel index */ 
j = jkc 7, nt; /•* dimension */ 
if ( j > nt2) j=nt-j ; 
w = j*dw; 

k = jkc/nt; /*ky dimension */ 
if (k > ny2) k=ny-k; 
ky = k*dky; 

dk = sqrt ( kx*kx + ky*ky ); 

if (dk > 0) v = v/dk; 
if (v > vtop | | dk == 0) { 
data0[jk] = 0; 
datal[jk] = 0; 

> 

/* if (k < 64) data0[jk] = 0; ♦/ 

} 


/•write the image back in fileO... •/ 

fsk = fseek(f ile[0] ,i*nf*sizeof (float) , SEEK. SET) ; 

cc = f write (fcdataO [0] , sizeof (float) ,nf ,file [0] ) ; 

if (cc ! = nf) IDL.MeasagedDL.M.NAMED. GENERIC, IDL.MSG.LONGJHP , 

"Write failure to file[0]."); 


/* . . .and f ilel */ 

fsk = fseek(filed] , (nx2-i-l)*nf*sizeof (float) .SEEK. SET) ; 

cc = fwrite(4datal[0] .sizeof (float), nf, filetl]); 

if (cc ! = nf) IDL.HessagedDL.M.NAMED.GEHERIC, IDL.MSG.LONC JMP , 

"Write failure to file[l]. M ); 
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/* Close the files */ 
for (i=0;i<NF;i++) { 
dum[0] = ftunit; 
unit. type = IDL.TYP.LONG; 
unit. flags = 0; 
unit. value. 1 = luns[i]; 
IDL.FileClose ( 1 , dum , NULL) ; 

> 


return tmpO; 

> 


Jtundef SWAP 
#undef NRANSI 



void f ilesetup( IDL.VPTR fnaaes ) 



{ 

int i, nfiles; 

IDL.FILE.STAT fstat; 

IDL.STRING *fn; 

IDL.VPTR dum[2] ; 

IDL.VARIABLE unit, name; 

/* Check properties »/ 

IDL.ENSURZ.STRING (f names) ; 

IDL.ENSURE. ARRAY (f names) ; 
if (fnames->value . arr->n_dim != 1) 

IDL.Message ( IDL.M.NAKED.CLYERIC , IDL.MSC.LONC JMP . 

"Temporary files array must be i-D. H ); 
nfiles = f names->value . arr->n elts; 
if (nfiles != NF) 

IDL.Message (IDL.M.NAMED.GEKERIC, IDL.MSC.LONGJMP, 

Wrong number of filenames passed."); 

/» printf ("\nf names type: 7.d\n" ,fnames->type) ; */ 

/ printf ( f names nelt = 7,d\n" , fnames->value . arr->n_elts) ; */ 

printf ("\nStorage Files:\n”); 

fn = &fnames->value.arr->data[0] ; 



/ 


/ 
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f or (i=0; i<nf iles ; i++) { 

/* printf ("filename 7d length = '/ t d\n" , i ,f n->slen) ; */ 

/* printf ("filename */.d type = ‘/.dXn", i ,fn->stype) ; */ 

printf ("File 7,d = */*s\n" t i >fn->s) ; 
fn++; 

> 

/* Open files */ 
printf ("\nOpening f iles. . . \n") ; 
fn = &fnames->value . arr->data[0] ; 
f or (i=0 ; i<nf iles ; i++) { 

/* get file unit • / 
dum[0] = frunit; 
unit. type = IDL.TYP.LOKG ; 
unit. flags = 0; 

IDL.FileGetUnit (1 ,dum) ; 
luns[i] = dum[0]->value.l; 

/* set up filename •/ 

/• printf ("Opening LUN Xd\n" l (int) luns[i]); •/ 

name. type = I DL.TYP. STRING; 
name. flags = IDL_V_CONST; 
name. value. str.s * fn->s; 
name. value. str.slen = sizeof (fn->s)-l ; 
name. value. str.stype = 0; 
dum[l] = Jcname ; 

/• NOTE: Extern. Dev. Guide function prototype incorrect: needs one more arg •/ 
IDL_File0pen(2, dum, (char •) 0, IDL.OPEN.R | IDL_0PEN_W, 

IDL.F_UNIX.F77 , 1,0); 

f n++ ; 

> 

printf ("Files successfully opened . \n" ) ; 

/•Check file status •/ 

printf ("\nEnsuring file status ... \n") ; 

for (i=0; i<nf iles; i++){ 

if (IDL.FileEnsureStatus (IDL.HSG.LONG JMP , (int) luns[i], 

IDL.EFS.OPEN | IDL.EFS.READ I IDL.EFS. WRITE) ) 
printfC'File '/.d is open for reading and writing. \n", 

(int) luns [i] ) ; 

} 

printf ("Status ensured. \n") ; 

printf ("\nGetting FILE pointers... \n"); 

for (i=0;i<nfile3;i++){ 
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IDL_FileStat((int) luns[i], fcfstat); 

/• printf ("Filename from FileStat() : */,s\n" .fstat .name) ; •/ 
file[i] = fstat. fptr; 

> 

printf ("FileStats complete. \n") ; 
return; 

> 
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